//
//  gl.c
//  gltest
//
//  Created by jiajiawei on 12-12-28.
//  Copyright (c) 2012年 jiajiawei. All rights reserved.
//

#include <stdio.h>
#include "gl.h"
#include "math.h"
#include <vector>
#include <algorithm>

#define ClipSegmentMaskAgree    0x0000
#define ClipSegmentMaskLeft     0x1000
#define ClipSegmentMaskTop      0x0100
#define ClipSegmentMaskRight    0x0010
#define ClipSegmentMaskBottom   0x0001


static unsigned int IMAGE_DATA[320 * 480];
static unsigned int TEXTURE_DATA[320 * 480];
static unsigned int DEPTHS_DATA[320 * 480];
static float MODEVIEW_M[16];
static float FRUSTUM_M[16];

unsigned int textureColor(float s, float t)
{
    if ( s < 0 )
    {
        s = -s;
    }
    if ( t < 0 )
    {
        t = -t;
    }
    s = roundf((s - int(s)) * 320);
    t = roundf((t - int(t)) * 480);

//    if ( s > 0.5 )
//    {
//        return 0xFF0000FF;
//    }
//    else
//    {
//        return 0xFF00FF00;
//    }
    return TEXTURE_DATA[(int)(t * 320 + s)];
}

void textureInit()
{
    UIImage* img = [UIImage imageWithContentsOfFile:[[NSBundle mainBundle] pathForResource:@"test" ofType:@"jpg"]];
    CGContextRef context = NULL;
    CGColorSpaceRef colorSpace = CGColorSpaceCreateDeviceRGB();
    context = CGBitmapContextCreate(TEXTURE_DATA, 320, 480, 8, 320 * 4, colorSpace, kCGImageAlphaPremultipliedLast | kCGBitmapByteOrder32Big);
    CGContextDrawImage(context, CGRectMake(0, 0, 320, 480), img.CGImage);
    CGColorSpaceRelease(colorSpace);
    CGContextRelease(context);
    return;
    for ( int i = 0; i < 320; i++ )
    {
        unsigned int color;
        for (int j = 0; j < 480; j++ )
        {
            if ( (j % 20) < 10 )
            {
                color = 0xFF0000FF;
            }
            else
            {
                color = 0xFFFF0000;
            }
            TEXTURE_DATA[j*320+i] = color;
        }

        
    }
}

unsigned int* image_data()
{
    return IMAGE_DATA;
}

unsigned int image_data_size()
{
    return sizeof(IMAGE_DATA);
}

float* modeviewM()
{
    return MODEVIEW_M;
}

float* frustumM()
{
    return FRUSTUM_M;
}



void setPixel(int x, int y, float z, unsigned int color)
{
    if ( x > 0 && x < 320 && y > 0 && y < 480 )
    {
        unsigned int cdepth = DEPTHS_DATA[y * 320 + x];
        unsigned int depth = (z + 1.0f) * UINT32_MAX / 2.0f;
        if ( depth < cdepth )
        {
            DEPTHS_DATA[y * 320 + x] = depth;
            IMAGE_DATA[y * 320 + x] = color;
        }
    }
}

unsigned int lineColor(int color1, int color2, float m)
{
    unsigned int color = 0;
    unsigned int r1 = (color1 >> 24) & 0xFF;
    unsigned int g1 = (color1 >> 16) & 0xFF;
    unsigned int b1 = (color1 >> 8 ) & 0xFF;
    unsigned int a1 = (color1 & 0xFF);
    unsigned int r2 = (color2 >> 24) & 0xFF;
    unsigned int g2 = (color2 >> 16) & 0xFF;
    unsigned int b2 = (color2 >> 8 ) & 0xFF;
    unsigned int a2 = (color2 & 0xFF);

    unsigned int  r = r1 * m + r2 * (1 - m);
    unsigned int  g = g1 * m + g2 * (1 - m);
    unsigned int  b = b1 * m + b2 * (1 - m);
    unsigned int  a = a1 * m + a2 * (1 - m);
    color += r << 24;
    color += g << 16;
    color += b << 8;
    color += a;
    return 0xFF0000FF;
}

inline int transWindowX(float x)
{
    return roundf(160 + x * 150);
}

inline int transWindowY(float y)
{
    return roundf(230 + y * 220);
}

inline float lerpf(float a, float b, float t)
{
    return a + (b - a) * t;
}

inline int lerpInt(int a, int b, float t)
{
    return roundf(a + t * (float)(b - a));
}

void genSLinePoins(GLPoint point0, GLPoint point1, std::vector<GLPoint>& points)
{
    int ax = roundf(point0.x);
    int ay = roundf(point0.y);
    int bx = roundf(point1.x);
    int by = roundf(point1.y);

    GLPoint pt;

    if ( ax == bx )
    {
        if ( ay > by )
        {
            std::swap(ay, by);
            std::swap(point0, point1);
        }
        for ( int i = ay; i <= by; i++ )
        {
            float kt = (float)(i - ay) / (float)(by - ay);
            pt.x = ax;
            pt.y = i;
            pt.z = lerpf(point0.z, point1.z, kt);
            pt.w = 1 / lerpf(1 / point0.w, 1 / point1.w, kt);
            pt.z = lerpf(point0.z / point0.w, point1.z / point1.w, kt) * pt.w;
            pt.s = lerpf(point0.s / point0.w, point1.s / point1.w, kt) * pt.w;
            pt.t = lerpf(point0.t / point0.w, point1.t / point1.w, kt) * pt.w;
//            pt.w = lerpf(point0.w, point1.w, kt);
//            kt = kt / lerpf(point1.w / point0.w, 1, kt);
//            pt.s = lerpf(point0.s, point1.s, kt);
//            pt.t = lerpf(point0.t, point1.t, kt);
//            float mu = lerpf(1.0f / point0.w, 1.0f / point1.w, kt);
//            pt.s = lerpf(point0.s / point0.w, point1.s / point1.w, kt) / mu;
//            pt.t = lerpf(point0.t / point0.w, point1.t / point1.w, kt) / mu;
            points.push_back(pt);
        }
        return;
    }

    if ( ay == by )
    {
        if ( ax > bx )
        {
            std::swap(ax, bx);
            std::swap(point0, point1);
        }
        for ( int i = ax; i <= bx; i++ )
        {
            float kt = (float)(i - ax) / (float)(bx - ax);

            pt.x = i;
            pt.y = ay;
            pt.z = lerpf(point0.z, point1.z, kt);
            pt.w = 1 / lerpf(1 / point0.w, 1 / point1.w, kt);
            pt.z = lerpf(point0.z / point0.w, point1.z / point1.w, kt) * pt.w;
            pt.s = lerpf(point0.s / point0.w, point1.s / point1.w, kt) * pt.w;
            pt.t = lerpf(point0.t / point0.w, point1.t / point1.w, kt) * pt.w;
//            pt.w = lerpf(point0.w, point1.w, kt);
//            kt = kt / lerpf(point1.w / point0.w, 1, kt);
//            pt.s = lerpf(point0.s, point1.s, kt);
//            pt.t = lerpf(point0.t, point1.t, kt);
//            float mu = lerpf(1.0f / point0.w, 1.0f / point1.w, kt);
//            pt.s = lerpf(point0.s / point0.w, point1.s / point1.w, kt) / mu;
//            pt.t = lerpf(point0.t / point0.w, point1.t / point1.w, kt) / mu;
            points.push_back(pt);
        }
        return;
    }

    float k = (float)(by - ay) / (float)(bx - ax);

    if ( k > 1 )
    {
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }
    else if ( k >= -1 && k < 0 )
    {
        ax = -ax;
        bx = -bx;
    }
    else if ( k < -1 )
    {
        ax = -ax;
        bx = -bx;
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }

    if ( ax > bx )
    {
        std::swap(ax, bx);
        std::swap(ay, by);
        std::swap(point0, point1);
    }

    int y = ay;
    int W = bx - ax;
    int H = by - ay;
    int F = 2 * H - W;
    for ( int x = ax; x <= bx; x++ )
    {
        if ( k > 1 )
        {
            pt = glPointMake(y, x, 0, 1);
        }
        else if ( k >= -1 && k < 0 )
        {
            pt = glPointMake(-x, y, 0, 1);
        }
        else if ( k < -1 )
        {
            pt = glPointMake(-y, x, 0, 1);
        }
        else
        {
            pt = glPointMake(x, y, 0, 1);
        }

        float kt = (float)(x - ax) / (float)(bx - ax);
        pt.z = lerpf(point0.z, point1.z, kt);
        pt.w = 1 / lerpf(1 / point0.w, 1 / point1.w, kt);
        pt.z = lerpf(point0.z / point0.w, point1.z / point1.w, kt) * pt.w;
        pt.s = lerpf(point0.s / point0.w, point1.s / point1.w, kt) * pt.w;
        pt.t = lerpf(point0.t / point0.w, point1.t / point1.w, kt) * pt.w;
//        pt.w = lerpf(point0.w, point1.w, kt);
//        kt = kt / lerpf(point1.w / point0.w, 1, kt);
//        pt.s = lerpf(point0.s, point1.s, kt);
//        pt.t = lerpf(point0.t, point1.t, kt);
//        float mu = lerpf(1.0f / point0.w, 1.0f / point1.w, kt);
//        pt.s = lerpf(point0.s / point0.w, point1.s / point1.w, kt) / mu;
//        pt.t = lerpf(point0.t / point0.w, point1.t / point1.w, kt) / mu;
        points.push_back(pt);


        if ( F < 0 )
        {
            F = F + 2 * H;
        }
        else
        {
            y = y + 1;
            F = F + 2 * (H - W);
        }
    }
    
    if ( k < 0 )
    {
        //std::reverse(vpoints.begin(), vpoints.end());
    }

}


void genLinePoints(std::vector<GLPoint>& vpoints)
{
    if ( vpoints.size() != 2 )
    {
        return;
    }

    int ax = roundf(vpoints.at(0).x);
    int ay = roundf(vpoints.at(0).y);
    int bx = roundf(vpoints.at(1).x);
    int by = roundf(vpoints.at(1).y);
    float at = vpoints.at(0).t;
    float as = vpoints.at(0).s;
    float bt = vpoints.at(1).t;
    float bs = vpoints.at(1).s;
    float t0 = at;
    float s0 = as;
    float dt = bt - at;
    float ds = bs - as;
    float t, s;
    float ad = vpoints.at(0).z;
    float bd = vpoints.at(1).z;
    float d0 = ad;
    float dd = bd - ad;
    float d;
    float aw = vpoints.at(0).w;
    float bw = vpoints.at(1).w;
    GLPoint pt;
    vpoints.clear();

    if ( ax == bx )
    {
        if ( ay > by )
        {
            int m;
            m = ay;
            ay = by;
            by = m;

            t0 = bt;
            s0 = bs;
            dt = -dt;
            ds = -ds;

            d0 = bd;
            dd = -dd;
        }
        for ( int i = ay; i <= by; i++ )
        {
            pt = glPointMake(ax, i, 0, 1);
            float kt = (float)(i - ay) / (float)(by - ay);
            s = (ds == 0 ? s0 : s0 + kt * ds);
            t = (dt == 0 ? t0 : t0 + kt * dt);
            d = (dd == 0 ? d0 : d0 + kt * dd);
            pt.t = t;
            pt.s = s;
            pt.z = d;
            vpoints.push_back(pt);
            
        }
        return;
    }

    if ( ay == by )
    {
        if ( ax > bx )
        {
            int m;
            m = ax;
            ax = bx;
            bx = m;

            t0 = bt;
            s0 = bs;
            dt = -dt;
            ds = -ds;

            d0 = bd;
            dd = -dd;
        }
        for ( int i = ax; i <= bx; i++ )
        {
            pt = glPointMake(i, ay, 0, 1);
            float kt = (float)(i - ax) / (float)(bx - ax);
            s = (ds == 0 ? s0 : s0 + kt * ds);
            t = (dt == 0 ? t0 : t0 + kt * dt);
            d = (dd == 0 ? d0 : d0 + kt * dd);
            pt.t = t;
            pt.s = s;
            pt.z = d;
            vpoints.push_back(pt);
        }
        return;
    }

    float k = (float)(by - ay) / (float)(bx - ax);

    if ( k > 1 )
    {
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }
    else if ( k >= -1 && k < 0 )
    {
        ax = -ax;
        bx = -bx;
    }
    else if ( k < -1 )
    {
        ax = -ax;
        bx = -bx;
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }

    if ( ax > bx )
    {
        int m;
        m = ax; ax = bx; bx = m;
        m = ay; ay = by; by = m;

        t0 = bt;
        s0 = bs;
        dt = -dt;
        ds = -ds;

        d0 = bd;
        dd = -dd;
    }

    int y = ay;
    int W = bx - ax;
    int H = by - ay;
    int F = 2 * H - W;
    for ( int x = ax; x <= bx; x++ )
    {
        if ( k > 1 )
        {
            pt = glPointMake(y, x, 0, 1);
        }
        else if ( k >= -1 && k < 0 )
        {
            pt = glPointMake(-x, y, 0, 1);
        }
        else if ( k < -1 )
        {
            pt = glPointMake(-y, x, 0, 1);
        }
        else
        {
            pt = glPointMake(x, y, 0, 1);
        }

        float kt = (float)(x - ax) / (float)(bx - ax);
        s = (ds == 0 ? s0 : s0 + kt * ds);
        t = (dt == 0 ? t0 : t0 + kt * dt);
        d = (dd == 0 ? d0 : d0 + kt * dd);
        pt.t = t;
        pt.s = s;
        pt.z = d;
        vpoints.push_back(pt);


        if ( F < 0 )
        {
            F = F + 2 * H;
        }
        else
        {
            y = y + 1;
            F = F + 2 * (H - W);
        }
    }

    if ( k < 0 )
    {
        //std::reverse(vpoints.begin(), vpoints.end());
    }
}


void drawLine(float x1, float y1, float x2, float y2, unsigned int color1, unsigned int color2)
{
    int ax = roundf(x1);
    int ay = roundf(y1);
    int bx = roundf(x2);
    int by = roundf(y2);
    if ( ax == bx )
    {
        if ( ay > by )
        {
            int m;
            m = ay;
            ay = by;
            by = m;

            m = color1;
            color1 = color2;
            color2 = m;
        }
        for ( int i = ay; i <= by; i++ )
        {
            setPixel(ax, i, 0, lineColor(color1, color2, (float)(i - ay) / (float)(by - ay)));
        }
        return;
    }

    if ( ay == by )
    {
        if ( ax > bx )
        {
            int m;
            m = ax;
            ax = bx;
            bx = m;
            m = color1;
            color1 = color2;
            color2 = m;
        }
        for ( int i = ax; i <= bx; i++ )
        {
            setPixel(i, ay, 0, lineColor(color1, color2, (float)(i - ax) / (float)(bx - ax)));
        }
        return;
    }

    float k = (float)(by - ay) / (float)(bx - ax);

    if ( k > 1 )
    {
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }
    else if ( k > -1 && k < 0 )
    {
        ax = -ax;
        bx = -bx;
    }
    else if ( k < -1 )
    {
        ax = -ax;
        bx = -bx;
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }

    if ( ax > bx )
    {
        int m;
        m = ax; ax = bx; bx = m;
        m = ay; ay = by; by = m;

        m = color1;
        color1 = color2;
        color2 = m;
    }

    int y = ay;
    int W = bx - ax;
    int H = by - ay;
    int F = 2 * H - W;
    for ( int x = ax; x <= bx; x++ )
    {
        unsigned int color = lineColor(color1, color2, (float)(x - ax) / (float)(bx - ax));
        if ( k > 1 )
        {
            setPixel(y, x, 0, color);
        }
        else if ( k > -1 && k < 0 )
        {
            setPixel(-x, y, 0, color);
        }
        else if ( k < -1 )
        {
            setPixel(-y, x, 0, color);
        }
        else
        {
            setPixel(x, y, 0, color);
        }
        if ( F < 0 )
        {
            F = F + 2 * H;
        }
        else
        {
            y = y + 1;
            F = F + 2 * (H - W);
        }
    }

}

typedef struct
{
    int xmin;
    int xmax;
    float sl;
    float sr;
    float tl;
    float tr;
    float dl;
    float dr;
    GLPoint xminPoint;
    GLPoint xmaxPoint;
}XDURATION;

int compareLessXDURATION(XDURATION& d1, XDURATION& d2)
{
    return d1.xmax < d2.xmax;
}

int compareLessSXDURATION(XDURATION& d1, XDURATION& d2)
{
    return d1.xmaxPoint.x < d2.xmaxPoint.x;
}




void drawTSTriagle(std::vector<GLPoint>& vpoints)
{
    typedef struct
    {
        int index;
        int count;
        int xmin;
        int xmax;
        int ymin;
        int ymax;
        int order;
    }LPT;


    if ( vpoints.size() < 3 )
    {
        return;
    }
    int lcount = vpoints.size();
    int xmin, xmax, ymin, ymax;
    
    std::vector<std::vector<GLPoint>> linepoints;
    std::vector<LPT> lpts;
    for ( int i = 0; i < lcount; i++ )
    {
        linepoints.push_back(std::vector<GLPoint>());
        lpts.push_back(LPT());
        std::vector<GLPoint>& tempPoints = linepoints.back();
        LPT& lpt = lpts.back();
        if ( i == (vpoints.size() - 1) )
        {
            genSLinePoins(vpoints.at(i), vpoints.at(0), tempPoints);
        }
        else
        {
            genSLinePoins(vpoints.at(i), vpoints.at(i + 1), tempPoints);
        }
        
        lpt.xmin = std::min(tempPoints.front().x, tempPoints.back().x);
        lpt.xmax = std::max(tempPoints.front().x, tempPoints.back().x);
        lpt.ymin = std::min(tempPoints.front().y, tempPoints.back().y);
        lpt.ymax = std::max(tempPoints.front().y, tempPoints.back().y);
        lpt.index = 0;
        lpt.count = tempPoints.size();
        
        if ( tempPoints.back().y > tempPoints.front().y )
        {
            lpt.order = 1;
        }
        else
        {
            lpt.order = -1;
        }

        if ( i == 0 )
        {
            xmin = lpt.xmin;
            xmax = lpt.xmax;
            ymin = lpt.ymin;
            ymax = lpt.ymax;
        }
        else
        {
            xmin = std::min(xmin, lpt.xmin);
            xmax = std::max(xmax, lpt.xmax);
            ymin = std::min(ymin, lpt.ymin);
            ymax = std::max(ymax, lpt.ymax);
        }
    }

    lcount = linepoints.size();
    for ( int y = ymin; y <= ymax; y++ )
    {
        std::vector<XDURATION> xpoints;
        XDURATION duration;
        for ( int i = 0; i < lcount; i++ )
        {
            LPT& lpt = lpts.at(i);
            if ( lpt.ymin == lpt.ymax )
            {
                continue;
            }
            if ( y < lpt.ymin || y > lpt.ymax || lpt.index >= lpt.count )
            {
                continue;
            }

            duration.xminPoint = duration.xmaxPoint = linepoints.at(i).at(lpt.index);
            lpt.index += 1;
            while (lpt.index < lpt.count && linepoints.at(i).at(lpt.index).y <= y)
            {
                duration.xmaxPoint = linepoints.at(i).at(lpt.index);
                lpt.index += 1;
            }

            if ( duration.xminPoint.x > duration.xmaxPoint.x )
            {
                std::swap(duration.xminPoint, duration.xmaxPoint);
            }

            //检测端点，去掉重复
            int nearIndex = -1;
            if ( y == lpt.ymin )
            {
                nearIndex = (lpt.order < 0 ? (i == (lcount - 1) ? 0 : i + 1) : (i == 0 ? lcount - 1 : i - 1));
            }
            if ( y == lpt.ymax )
            {
                nearIndex = (lpt.order > 0 ? (i == (lcount - 1) ? 0 : i + 1) : (i == 0 ? lcount - 1 : i - 1));
            }
            if ( nearIndex != -1 )
            {
                LPT& nearlpt = lpts.at(nearIndex);
                if ( nearlpt.ymax != nearlpt.ymin && nearIndex > i )
                {
                    continue;
                }
            }
            xpoints.push_back(duration);
        }
        std::sort(xpoints.begin(), xpoints.end(), compareLessSXDURATION);
        for ( int index = 0; index < (xpoints.size()/2); index++)
        {
            XDURATION& d0 = xpoints.at(index);
            XDURATION& d1 = xpoints.at(index + 1);
            GLPoint xminPoint, xmaxPoint;

            xminPoint = (d0.xminPoint.x < d1.xminPoint.x ? d0.xminPoint : d1.xminPoint);
            xmaxPoint = (d0.xmaxPoint.x > d1.xminPoint.x ? d0.xmaxPoint : d1.xmaxPoint);

            for ( int x = xminPoint.x; x <= xmaxPoint.x; x++)
            {
                float kt = (float)(x - xminPoint.x) / (float)(xmaxPoint.x - xminPoint.x);
                float s, t, z;

                z = lerpf(xminPoint.z, xmaxPoint.z, kt);
                float w = 1 / lerpf(1 / xminPoint.w, 1 / xmaxPoint.w, kt);
                s = lerpf(xminPoint.s / xminPoint.w, xmaxPoint.s / xmaxPoint.w, kt) * w;
                t = lerpf(xminPoint.t / xminPoint.w, xmaxPoint.t / xmaxPoint.w, kt) * w;
                z = lerpf(xminPoint.z / xminPoint.w, xmaxPoint.z / xmaxPoint.w, kt) * w;
                setPixel(x, y, z / w, textureColor(s, t));
            }
        }
    }
}


void drawSTriagle(std::vector<GLPoint>& vpoints)
{
    typedef struct
    {
        int index;
        int count;
        int xmin;
        int xmax;
        int ymin;
        int ymax;
        int order;
    }LPT;

    
    if ( vpoints.size() < 3 )
    {
        return;
    }
    int lcount = vpoints.size();
    int xmin, xmax, ymin, ymax;
    std::vector<GLPoint> tpoints;
    std::vector<std::vector<GLPoint>> linepoints;
    std::vector<LPT> lpts;
    for ( int i = 0; i < lcount; i++ )
    {
        tpoints.clear();
        if ( i == (vpoints.size() - 1) )
        {
            tpoints.push_back(vpoints.at(i));
            tpoints.push_back(vpoints.at(0));
        }
        else
        {
            tpoints.push_back(vpoints.at(i));
            tpoints.push_back(vpoints.at(i+1));
        }
        if ( tpoints.front().y == tpoints.back().y )
        {
            //continue;
        }
        linepoints.push_back(tpoints);
        genLinePoints(linepoints.back());

        std::vector<GLPoint>& tempPoints = linepoints.back();
        LPT lpt;
        lpt.xmin = std::min(tempPoints.front().x, tempPoints.back().x);
        lpt.xmax = std::max(tempPoints.front().x, tempPoints.back().x);
        lpt.ymin = std::min(tempPoints.front().y, tempPoints.back().y);
        lpt.ymax = std::max(tempPoints.front().y, tempPoints.back().y);
        lpt.index = 0;
        lpt.count = tempPoints.size();
        if ( tpoints.back().y > tpoints.front().y )
        {
            lpt.order = 1;
        }
        else
        {
            lpt.order = -1;
        }
        lpts.push_back(lpt);

        if ( i == 0 )
        {
            xmin = lpt.xmin;
            xmax = lpt.xmax;
            ymin = lpt.ymin;
            ymax = lpt.ymax;
        }
        else
        {
            xmin = std::min(xmin, lpt.xmin);
            xmax = std::max(xmax, lpt.xmax);
            ymin = std::min(ymin, lpt.ymin);
            ymax = std::max(ymax, lpt.ymax);
        }
    }

    lcount = linepoints.size();
    for ( int y = ymin; y <= ymax; y++ )
    {
        std::vector<XDURATION> xpoints;
        for ( int i = 0; i < lcount; i++ )
        {
            LPT& lpt = lpts.at(i);
            if ( lpt.ymin == lpt.ymax )
            {
                continue;
            }
            if ( y >= lpt.ymin && y <= lpt.ymax && lpt.index < lpt.count )
            {
                int dxmin, dxmax;
                float sl, sr, tl, tr;
                float dl, dr;
                dxmax = dxmin = linepoints.at(i).at(lpt.index).x;
                sl = sr = linepoints.at(i).at(lpt.index).s;
                tl = tr = linepoints.at(i).at(lpt.index).t;
                dl = dr = linepoints.at(i).at(lpt.index).z;
                lpt.index += 1;
                while (lpt.index < lpt.count && linepoints.at(i).at(lpt.index).y <= y)
                {
                    dxmax = linepoints.at(i).at(lpt.index).x;
                    sr = linepoints.at(i).at(lpt.index).s;
                    tr = linepoints.at(i).at(lpt.index).t;
                    dr = linepoints.at(i).at(lpt.index).z;
                    lpt.index += 1;
                }
                
                if ( dxmin > dxmax )
                {
                    std::swap(dxmin, dxmax);
                    std::swap(tl, tr);
                    std::swap(sl, sr);
                    std::swap(dl, dr);
                }

                //检测端点，去掉重复
                if ( y == lpt.ymin )
                {
                    int nearIndex = -1;
                    if ( lpt.order < 0 )
                    {
                        if ( i == (lcount - 1) )
                        {
                            nearIndex = 0;
                        }
                        else
                        {
                            nearIndex = i + 1;
                        }
                    }
                    else
                    {
                        if ( i == 0 )
                        {
                            nearIndex = lcount - 1;
                        }
                        else
                        {
                            nearIndex = i - 1;
                        }
                    }
                    
                    LPT& nearlpt = lpts.at(nearIndex);
                    if ( nearlpt.ymin < y && nearlpt.ymin != nearlpt.ymax )
                    {
                        continue;
                    }
                }

                XDURATION xduration;
                xduration.xmin = dxmin;
                xduration.xmax = dxmax;
                xduration.sl = sl;
                xduration.sr = sr;
                xduration.tl = tl;
                xduration.tr = tr;
                xduration.dl = dl;
                xduration.dr = dr;
                xpoints.push_back(xduration);
            }
        }
        std::sort(xpoints.begin(), xpoints.end(), compareLessXDURATION);
        if ( (xpoints.size() % 2) == 1 )
        {
            NSLog(@"count=%ld", xpoints.size());
        }
        for ( int index = 0; index < (xpoints.size()/2); index++)
        {
            XDURATION& d0 = xpoints.at(index);
            XDURATION& d1 = xpoints.at(index + 1);
            int dxmin, dxmax;
            float s0, s1, t0, t1, s, t;
            float z0, z1, z;
            if ( d0.xmin < d1.xmin )
            {
                dxmin = d0.xmin;
                s0 = d0.sl;
                t0 = d0.tl;
                z0 = d0.dl;
            }
            else
            {
                dxmin = d1.xmin;
                s0 = d1.sl;
                t0 = d1.tl;
                z0 = d1.dl;
            }

            if ( d0.xmax > d1.xmax )
            {
                dxmax = d0.xmax;
                s1 = d0.sr;
                t1 = d0.tr;
                z1 = d0.dr;
            }
            else
            {
                dxmax = d1.xmax;
                s1 = d1.sr;
                t1 = d1.tr;
                z1 = d1.dr;
            }
            
            
            for ( int x = dxmin; x <= dxmax; x++)
            {
                float kt = (float)(x - dxmin) / (float)(dxmax - dxmin);
                s = (s1 == s0 ? s0 : s0 + kt * (s1 - s0));
                t = (t1 == t0 ? t0 : t0 + kt * (t1 - t0));
                z = (z1 == z0 ? z0 : z0 + kt * (z1 - z0));
                setPixel(x, y, z, textureColor(s, t));
            }
        }
    }
}

void drawTriagle(float x1, float y1, float x2, float y2, float x3, float y3, unsigned int color1, unsigned int color2, unsigned int color3)
{
    


    int ax = roundf(x1);
    int ay = roundf(y1);
    int bx = roundf(x2);
    int by = roundf(y2);
    int cx = roundf(x3);
    int cy = roundf(y3);
    if ( ax == bx )
    {
        if ( ay > by )
        {
            int m;
            m = ay;
            ay = by;
            by = m;
            m = color1;
            color1 = color2;
            color2 = m;
        }
        for ( int i = ay; i <= by; i++ )
        {
            unsigned int color = lineColor(color1, color2, (float)(i - ay) / (float)(by - ay));
            drawLine(ax, i, cx, cy, color, color3);
            //setPixel(ax, i, color);
        }
        return;
    }

    if ( ay == by )
    {
        if ( ax > bx )
        {
            int m;
            m = ax;
            ax = ay;
            ay = m;
            m = color1;
            color1 = color2;
            color2 = m;
        }
        for ( int i = ax; i <= bx; i++ )
        {
            unsigned int color = lineColor(color1, color2, (float)(i - ax) / (float)(bx - ax));
            drawLine(i, ay, cx, cy, color, color3);
            //setPixel(i, ay, color);
        }
        return;
    }

    float k = (float)(by - ay) / (float)(bx - ax);

    if ( k > 1 )
    {
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }
    else if ( k > -1 && k < 0 )
    {
        ax = -ax;
        bx = -bx;
    }
    else if ( k < -1 )
    {
        ax = -ax;
        bx = -bx;
        int m;
        m = ax; ax = ay; ay = m;
        m = bx; bx = by; by = m;
    }

    if ( ax > bx )
    {
        int m;
        m = ax; ax = bx; bx = m;
        m = ay; ay = by; by = m;

        m = color1;
        color1 = color2;
        color2 = m;
    }

    int y = ay;
    int W = bx - ax;
    int H = by - ay;
    int F = 2 * H - W;
    for ( int x = ax; x <= bx; x++ )
    {
        unsigned int color = lineColor(color1, color2, (float)(x - ax) / (float)(bx - ax));
        if ( k > 1 )
        {
            drawLine(y, x, cx, cy, color, color3);
            //setPixel(y, x, color);
        }
        else if ( k > -1 && k < 0 )
        {
            drawLine(-x, y, cx, cy, color, color3);
            //setPixel(-x, y, color);
        }
        else if ( k < -1 )
        {
            drawLine(-y, x, cx, cy, color, color3);
            //setPixel(-y, x, color);
        }
        else
        {
            drawLine(x, y, cx, cy, color, color3);
            //setPixel(x, y, color);
        }
        if ( F < 0 )
        {
            F = F + 2 * H;
        }
        else
        {
            y = y + 1;
            F = F + 2 * (H - W);
        }
    }
}



GLPoint glPointMake(float x, float y, float z, float w)
{
    GLPoint point;
    point.x = x;
    point.y = y;
    point.z = z;
    point.w = w;
    return point;
}

void glPointSet(GLPoint* point, float x, float y, float z, float w)
{
    point->x = x;
    point->y = y;
    point->z = z;
    point->w = w;
}

void glPointNormal(GLPoint* point)
{
    point->x = point->x / point->w;
    point->y = point->y / point->w;
    point->z = point->z / point->w;
    point->w = 1.0f;
}

GLRect GLRectMake(float left, float top, float right, float bottom)
{
    GLRect rect;
    rect.left = left;
    rect.top = top;
    rect.right = right;
    rect.bottom = bottom;
    return rect;
}

GLVector glVectorMake(float x, float y, float z)
{
    GLVector vector;
    vector.x = x;
    vector.y = y;
    vector.z = z;
    vector.n = 0;
    return vector;
}

void glVectorSet(GLVector* v, float x, float y, float z)
{
    v->x = x;
    v->y = y;
    v->z = z;
    v->n = 0;
}

void glVectorSub(GLVector* v1, GLVector* v2, GLVector* v)
{
    v->x = v1->x - v2->x;
    v->y = v1->y - v2->y;
    v->z = v1->z - v2->z;
    v->n = 0;
}

void glNormalVector(GLVector* v)
{
    float s = sqrt(v->x * v->x + v->y * v->y + v->z * v->z);
    if ( s == 0.0f )
    {
        return;
    }
    v->x = v->x / s;
    v->y = v->y / s;
    v->z = v->z / s;
    v->n = 0;
}

float glVectorMul(GLVector* v1, GLVector* v2)
{
    return v1->x * v2->x + v1->y * v2->y + v1->z * v2->z;
}

void glVectorCross(GLVector* v1, GLVector* v2, GLVector* v)
{
    v->x = v1->y * v2->z - v1->z * v2->y;
    v->y = -v1->x * v2->z + v1->z * v2->x;
    v->z = v1->x * v2->y - v1->y * v2->x;
    v->n = 0;
}

int clipSegmentMask(GLPoint* point, GLRect* window)
{
    int mask = 0;
    if ( point->x < window->left )
    {
        mask = mask | ClipSegmentMaskLeft;
    }
    if ( point->y < window->top )
    {
        mask = mask | ClipSegmentMaskTop;
    }
    if ( point->x > window->right )
    {
        mask = mask | ClipSegmentMaskRight;
    }
    if ( point->y > window->bottom )
    {
        mask = mask | ClipSegmentMaskBottom;
    }
    return mask;
}

void clipPoint(GLPoint* point, GLPoint point0, int* mask, int delx, int dely, GLRect* window)
{
    if ( (*mask & ClipSegmentMaskLeft) == ClipSegmentMaskLeft )
    {
        point->x = window->left;
        point->y = point0.y + (window->left - point0.x) * dely / delx;
    }
    else if ( (*mask & ClipSegmentMaskTop) == ClipSegmentMaskTop )
    {
        point->y = window->top;
        point->x = point0.x + (window->top - point0.y) * delx / dely;
    }
    else if ( (*mask & ClipSegmentMaskRight) == ClipSegmentMaskRight )
    {
        point->x = window->right;
        point->y = point0.y + (window->right - point0.x) * dely / delx;
    }
    else if ( (*mask & ClipSegmentMaskBottom) == ClipSegmentMaskBottom )
    {
        point->y = window->bottom;
        point->x = point0.x + (window->bottom - point0.y) * delx / dely;
    }
    *mask = clipSegmentMask(point, window);
}

int clipSegment(GLPoint* point1, GLPoint* point2, GLRect* window)
{
    int mask1 = clipSegmentMask(point1, window);
    int mask2 = clipSegmentMask(point2, window);

    while ( 1 )
    {
        if ( mask1 == ClipSegmentMaskAgree && mask2 == ClipSegmentMaskAgree )
        {
            return 1;
        }

        if ( (mask1 & mask2) != 0 )
        {
            return 0;
        }

        float delx = point2->x - point1->x;
        float dely = point2->y - point1->y;

        if ( mask1 > 0 )
        {
            clipPoint(point1, *point1, &mask1, delx, dely, window);
        }
        else
        {
            clipPoint(point2, *point1, &mask2, delx, dely, window);
        }
    }
}

int computeBC(GLPoint* point, float pointBC[6])
{
    int outcode = 0;
    pointBC[0] = point->w + point->x;
    pointBC[1] = point->w - point->x;
    pointBC[2] = point->w + point->y;
    pointBC[3] = point->w - point->y;
    pointBC[4] = point->w + point->z;
    pointBC[5] = point->w - point->z;

    for (int i = 0; i < 6; i++ )
    {
        if ( pointBC[i] < 0 )
        {
            outcode = outcode | (0x1 << i);
        }
    }
    return outcode;
}

int clipEdge(GLPoint* point1, GLPoint* point2)
{
    float tIn = 0, tOut = 1.0f, tHit;
    float p1BC[6], p2BC[6];
    int p1Outcode, p2Outcode;
    p1Outcode = computeBC(point1, p1BC);
    p2Outcode = computeBC(point2, p2BC);

    if ( (p1Outcode & p2Outcode) != 0 )
    {
        return 0;
    }

    if ( (p1Outcode | p2Outcode) == 0 )
    {
        return 1;
    }

    for ( int i = 0; i < 6; i++ )
    {
        if ( p2BC[i] < 0 )
        {
            tHit = p1BC[i] / (p1BC[i] - p2BC[i]);
            tOut = fminf(tOut, tHit);
        }
        else if( p1BC[i] < 0 )
        {
            tHit = tHit = p1BC[i] / (p1BC[i] - p2BC[i]);
            tIn = fmaxf(tIn, tHit);
        }

        if ( tIn > tOut )
        {
            return 0;
        }
    }

    GLPoint tmp = *point1;
    if ( p1Outcode != 0 )
    {
        tmp.x = point1->x + tIn * (point2->x - point1->x);
        tmp.y = point1->y + tIn * (point2->y - point1->y);
        tmp.z = point1->z + tIn * (point2->z - point1->z);
        tmp.w = point1->w + tIn * (point2->w - point1->w);
    }

    if ( p2Outcode != 0 )
    {
        point2->x = point1->x + tOut * (point2->x - point1->x);
        point2->y = point1->y + tOut * (point2->y - point1->y);
        point2->z = point1->z + tOut * (point2->z - point1->z);
        point2->w = point1->w + tOut * (point2->w - point1->w);
    }

    glPointSet(point1, tmp.x, tmp.y, tmp.z, tmp.w);
    return 1;
}

void insertPoint(GLPoint* points, int* count, GLPoint point, int index)
{
    for ( int i = *count; i > index; i-- )
    {
        points[i] = points[i-1];
    }
    points[index] = point;
    *count += 1;
}

void removePoint(GLPoint* points, int* count, int index)
{
    for ( int i = index; i < *count - 1; i++ )
    {
        points[i] = points[i+1];
    }
    *count -= 1;
}


#define FLAGLEFT(point) (point.w + point.x)
#define FLAGRIGHT(point) (point.w - point.x)
#define FLAGBACK(point) (point.w + point.y)
#define FLAGFRONT(point) (point.w - point.y)
#define FLAGBOTTOM(point) (point.w + point.z)
#define FLAGTOP(point) (point.w - point.z)

#define TLEFT(point1, point2) (-(point1.w + point1.x) / ((point2.w - point1.w) + (point2.x - point1.x)))
#define TRIGHT(point1, point2) (-(point1.w - point1.x) / ((point2.w - point1.w) - (point2.x - point1.x)))
#define TBACK(point1, point2) (-(point1.w + point1.y) / ((point2.w - point1.w) + (point2.y - point1.y)))
#define TFRONT(point1, point2) (-(point1.w - point1.y) / ((point2.w - point1.w) - (point2.y - point1.y)))
#define TBOTTOM(point1, point2) (-(point1.w + point1.z) / ((point2.w - point1.z) + (point2.x - point1.z)))
#define TTOP(point1, point2) (-(point1.w - point1.z) / ((point2.w - point1.z) - (point2.x - point1.z)))

#define TPOINT(point1, point2, point, t)\
{\
    float f = point2.w * t / ( point1.w * (1-t) + point2.w * t);\
    point.x = point1.x + t * (point2.x - point1.x);\
    point.y = point1.y + t * (point2.y - point1.y);\
    point.z = point1.z + f * (point2.z - point1.z) * (point1.w + t * (point2.w - point1.w));\
    point.w = point1.w + t * (point2.w - point1.w);\
    point.s = point1.s + t * (point2.s - point1.s);\
    point.t = point1.t + t * (point2.t - point1.t);\
}

//point.z = -(point2.w * t) / (point1.w * (1-t) + point2.w * t);\

int clipTriagle(GLPoint* points, int* count)
{
    int index, len;

#define VFLAG FLAGLEFT
#define VT TLEFT
    index = 0;
    len = *count;
    for ( int i = 0; i < len; i++ )
    {
        GLPoint point = points[index], lpoint, rpoint;
        float flag = VFLAG(point);
        if ( flag < 0 )
        {
            lpoint = (index > 0 ? points[index - 1] : points[*count - 1]);
            rpoint = (index < (*count - 1) ? points[index + 1] : points[0]);
            float lflag = VFLAG(lpoint);
            float rflag = VFLAG(rpoint);
            if ( lflag <= 0 && rflag <= 0 )
            {
                removePoint(points, count, index);
            }
            else if ( lflag > 0 || rflag > 0 )
            {
                removePoint(points, count, index);
                GLPoint tpoint;
                float t;
                if ( lflag > 0 )
                {
                    t = VT(point, lpoint);
                    TPOINT(point, lpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
                if ( rflag > 0 )
                {
                    t = VT(point, rpoint);
                    TPOINT(point, rpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
            }
        }
        else
        {
            index += 1;
        }
    }
#undef VFLAG
#undef VT


//right
#define VFLAG FLAGRIGHT
#define VT TRIGHT
    index = 0;
    len = *count;
    for ( int i = 0; i < len; i++ )
    {
        GLPoint point = points[index], lpoint, rpoint;
        float flag = VFLAG(point);
        if ( flag < 0 )
        {
            lpoint = (index > 0 ? points[index - 1] : points[*count - 1]);
            rpoint = (index < (*count - 1) ? points[index + 1] : points[0]);
            float lflag = VFLAG(lpoint);
            float rflag = VFLAG(rpoint);
            if ( lflag <= 0 && rflag <= 0 )
            {
                removePoint(points, count, index);
            }
            else if ( lflag > 0 || rflag > 0 )
            {
                removePoint(points, count, index);
                GLPoint tpoint;
                float t;
                if ( lflag > 0 )
                {
                    t = VT(point, lpoint);
                    TPOINT(point, lpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
                if ( rflag > 0 )
                {
                    t = VT(point, rpoint);
                    TPOINT(point, rpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
            }
        }
        else
        {
            index += 1;
        }
    }
#undef VFLAG
#undef VT


//back
#define VFLAG FLAGBACK
#define VT TBACK
    index = 0;
    len = *count;
    for ( int i = 0; i < len; i++ )
    {
        GLPoint point = points[index], lpoint, rpoint;
        float flag = VFLAG(point);
        if ( flag < 0 )
        {
            lpoint = (index > 0 ? points[index - 1] : points[*count - 1]);
            rpoint = (index < (*count - 1) ? points[index + 1] : points[0]);
            float lflag = VFLAG(lpoint);
            float rflag = VFLAG(rpoint);
            if ( lflag <= 0 && rflag <= 0 )
            {
                removePoint(points, count, index);
            }
            else if ( lflag > 0 || rflag > 0 )
            {
                removePoint(points, count, index);
                GLPoint tpoint;
                float t;
                if ( lflag > 0 )
                {
                    t = VT(point, lpoint);
                    TPOINT(point, lpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
                if ( rflag > 0 )
                {
                    t = VT(point, rpoint);
                    TPOINT(point, rpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
            }
        }
        else
        {
            index += 1;
        }
    }
#undef VFLAG
#undef VT

//front
#define VFLAG FLAGFRONT
#define VT TFRONT
    index = 0;
    len = *count;
    for ( int i = 0; i < len; i++ )
    {
        GLPoint point = points[index], lpoint, rpoint;
        float flag = VFLAG(point);
        if ( flag < 0 )
        {
            lpoint = (index > 0 ? points[index - 1] : points[*count - 1]);
            rpoint = (index < (*count - 1) ? points[index + 1] : points[0]);
            float lflag = VFLAG(lpoint);
            float rflag = VFLAG(rpoint);
            if ( lflag <= 0 && rflag <= 0 )
            {
                removePoint(points, count, index);
            }
            else if ( lflag > 0 || rflag > 0 )
            {
                removePoint(points, count, index);
                GLPoint tpoint;
                float t;
                if ( lflag > 0 )
                {
                    t = VT(point, lpoint);
                    TPOINT(point, lpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
                if ( rflag > 0 )
                {
                    t = VT(point, rpoint);
                    TPOINT(point, rpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
            }
        }
        else
        {
            index += 1;
        }
    }
#undef VFLAG
#undef VT


//bottom
#define VFLAG FLAGBOTTOM
#define VT TBOTTOM
    index = 0;
    len = *count;
    for ( int i = 0; i < len; i++ )
    {
        GLPoint point = points[index], lpoint, rpoint;
        float flag = VFLAG(point);
        if ( flag < 0 )
        {
            lpoint = (index > 0 ? points[index - 1] : points[*count - 1]);
            rpoint = (index < (*count - 1) ? points[index + 1] : points[0]);
            float lflag = VFLAG(lpoint);
            float rflag = VFLAG(rpoint);
            if ( lflag <= 0 && rflag <= 0 )
            {
                removePoint(points, count, index);
            }
            else if ( lflag > 0 || rflag > 0 )
            {
                removePoint(points, count, index);
                GLPoint tpoint;
                float t;
                if ( lflag > 0 )
                {
                    t = VT(point, lpoint);
                    TPOINT(point, lpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
                if ( rflag > 0 )
                {
                    t = VT(point, rpoint);
                    TPOINT(point, rpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
            }
        }
        else
        {
            index += 1;
        }
    }
#undef VFLAG
#undef VT

//top
#define VFLAG FLAGTOP
#define VT TTOP
    index = 0;
    len = *count;
    for ( int i = 0; i < len; i++ )
    {
        GLPoint point = points[index], lpoint, rpoint;
        float flag = VFLAG(point);
        if ( flag < 0 )
        {
            lpoint = (index > 0 ? points[index - 1] : points[*count - 1]);
            rpoint = (index < (*count - 1) ? points[index + 1] : points[0]);
            float lflag = VFLAG(lpoint);
            float rflag = VFLAG(rpoint);
            if ( lflag <= 0 && rflag <= 0 )
            {
                removePoint(points, count, index);
            }
            else
            {
                removePoint(points, count, index);
                GLPoint tpoint;
                float t;
                if ( lflag > 0 )
                {
                    t = VT(point, lpoint);
                    TPOINT(point, lpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
                if ( rflag > 0 )
                {
                    t = VT(point, rpoint);
                    TPOINT(point, rpoint, tpoint, t);
                    insertPoint(points, count, tpoint, index);
                    index += 1;
                }
            }
        }
        else
        {
            index += 1;
        }
    }
#undef VFLAG
#undef VT

    return 1;
}


void glLoadIdentifier()
{
    MODEVIEW_M[0] = 1.0f;
    MODEVIEW_M[1] = 0.0f;
    MODEVIEW_M[2] = 0.0f;
    MODEVIEW_M[3] = 0.0f;
    
    MODEVIEW_M[4] = 0.0f;
    MODEVIEW_M[5] = 1.0f;
    MODEVIEW_M[6] = 0.0f;
    MODEVIEW_M[7] = 0.0f;
    
    MODEVIEW_M[8] = 0.0f;
    MODEVIEW_M[9] = 0.0f;
    MODEVIEW_M[10] = 1.0f;
    MODEVIEW_M[11] = 0.0f;

    MODEVIEW_M[12] = 0.0f;
    MODEVIEW_M[13] = 0.0f;
    MODEVIEW_M[14] = 0.0f;
    MODEVIEW_M[15] = 1.0f;
}

void glLoadMatrix(float m[16])
{
    for ( int i = 0; i < 16; i++ )
    {
        MODEVIEW_M[i] = m[i];
    }
}

void glMatrixMul(float m1[16], float m2[16], float m[16])
{
    float TM[16];
    for (int i = 0; i < 4; i++)
    {
        for (int j = 0; j < 4; j++)
        {
            TM[i + 4 * j] = m1[i] * m2[4*j] + m1[i+4] * m2[4*j+1] + m1[i+8] * m2[4*j+2] + m1[i+12] * m2[4*j+3];
        }
    }
    for ( int i = 0; i < 16; i++ )
    {
        m[i] = TM[i];
    }
}


void glMulMatrix(float m[16])
{
    float TM[16];
    for (int i = 0; i < 4; i++)
    {
        for (int j = 0; j < 4; j++)
        {
            TM[i + 4 * j] = MODEVIEW_M[i] * m[4*j] + MODEVIEW_M[i+4] * m[4*j+1] + MODEVIEW_M[i+8] * m[4*j+2] + MODEVIEW_M[i+12] * m[4*j+3];
        }
    }
    for ( int i = 0; i < 16; i++ )
    {
        MODEVIEW_M[i] = TM[i];
    }
}

void glTranslate(float x, float y, float z)
{
    float TM[16] =
    {
        1, 0, 0, 0,
        0, 1, 0, 0,
        0, 0, 1, 0,
        x, y, z, 1
    };
    glMulMatrix(TM);
}

void glScale(float x, float y, float z)
{
    float TM[16] =
    {
        x, 0, 0, 0,
        0, y, 0, 0,
        0, 0, z, 0,
        0, 0, 0, 1
    };
    glMulMatrix(TM);
}

void glRotate(float angle, float x, float y, float z)
{
    float tangle = GL_PI * angle / 180.0f;
    float sint = sinf(tangle);
    float cst = cosf(tangle);
    float TM[16] =
    {
        (float)(cst + (1 - cst) * pow(x, 2)), (1 - cst) * x * y + sint * z, (1 - cst) * x * z - sint * y, 0,
        (1 - cst) * x * y - sint * z, (float)(cst + (1 - cst) * pow(y, 2)), (1 - cst) * z * y + sint * x, 0,
        (1 - cst) * x * z + sint * y, (1 - cst) * z * y - sint * x, (float)(cst + (1 - cst) * pow(z, 2)), 0,
        0, 0, 0, 1
    };
    glMulMatrix(TM);
}

void glTransformPoint(float matrix[16], GLPoint* point)
{
    float px, py, pz, pg;
    px = matrix[0] * point->x + matrix[4] * point->y + matrix[8] * point->z + matrix[12] * point->w;
    py = matrix[1] * point->x + matrix[5] * point->y + matrix[9] * point->z + matrix[13] * point->w;
    pz = matrix[2] * point->x + matrix[6] * point->y + matrix[10] * point->z + matrix[14] * point->w;
    pg = matrix[3] * point->x + matrix[7] * point->y + matrix[11] * point->z + matrix[15] * point->w;
    point->x = px;
    point->y = py;
    point->z = pz;
    point->w = pg;
}

void gluLookAt(float eyex, float eyey, float eyez, float lookx, float looky, float lookz, float upx, float upy, float upz)
{
    GLVector eye = glVectorMake(eyex, eyey, eyez);
    GLVector look = glVectorMake(lookx, looky, lookz);
    GLVector up = glVectorMake(upx, upy, upz);
    
    GLVector n, u, v;
    glVectorSub(&eye, &look, &n);
    glVectorCross(&up, &n, &u);
    glVectorCross(&n, &u, &v);
    glNormalVector(&n);
    glNormalVector(&u);
    glNormalVector(&v);
    float dx, dy, dz;
    dx = -glVectorMul(&eye, &u);
    dy = -glVectorMul(&eye, &v);
    dz = -glVectorMul(&eye, &n);
    float V[16] =
    {
        u.x,    v.x,    n.x,    0,
        u.y,    v.y,    n.y,    0,
        u.z,    v.z,    n.z,    0,
        dx,     dy,     dz,     1
    };
    glMulMatrix(V);
}

void glCameraSetCT(GLCamera* camera)
{
    GLVector veye = glVectorMake(camera->eye.x, camera->eye.y, camera->eye.z);
    camera->CT[0] = camera->u.x;
    camera->CT[4] = camera->u.y;
    camera->CT[8] = camera->u.z;
    camera->CT[12] = -glVectorMul(&veye, &camera->u);
    camera->CT[1] = camera->v.x;
    camera->CT[5] = camera->v.y;
    camera->CT[9] = camera->v.z;
    camera->CT[13] = -glVectorMul(&veye, &camera->v);
    camera->CT[2] = camera->n.x;
    camera->CT[6] = camera->n.y;
    camera->CT[10] = camera->n.z;
    camera->CT[14] = -glVectorMul(&veye, &camera->n);
    camera->CT[3] = 0;
    camera->CT[7] = 0;
    camera->CT[11] = 0;
    camera->CT[15] = 1;
}

void glCameraSet(GLCamera* camera, float eyex, float eyey, float eyez, float lookx, float looky, float lookz, float upx, float upy, float upz)
{
    glPointSet(&camera->eye, eyex, eyey, eyez, 1);
    glPointSet(&camera->look, lookx, looky, lookz, 1);
    glVectorSet(&camera->up, upx, upy, upz);

    GLVector veye = glVectorMake(camera->eye.x, camera->eye.y, camera->eye.z);
    GLVector vlook = glVectorMake(camera->look.x, camera->look.y, camera->look.z);
    glVectorSub(&veye, &vlook, &camera->n);
    glVectorCross(&camera->up, &camera->n, &camera->u);
    glVectorCross(&camera->n, &camera->u, &camera->v);
    glNormalVector(&camera->u);
    glNormalVector(&camera->v);
    glNormalVector(&camera->n);
    glCameraSetCT(camera);
}

void glCameraSlide(GLCamera* camera, float du, float dv, float dn)
{
    float chgx = du * camera->u.x + dv * camera->v.x + dn * camera->n.x;
    float chgy = du * camera->u.y + dv * camera->v.y + dn * camera->n.y;
    float chgz = du * camera->u.z + dv * camera->v.z + dn * camera->n.z;

    glPointSet(&camera->eye, camera->eye.x + chgx, camera->eye.y + chgy, camera->eye.z + chgz, 1);
    glPointSet(&camera->look, camera->look.x + chgx, camera->look.y + chgy, camera->look.z + chgz, 1);
    glCameraSetCT(camera);
}

void glCameraRoll(GLCamera* camera, float angle)
{
    float tangle = GL_PI * angle / 180;
    float cs = cosf(tangle);
    float sn = sinf(tangle);
    GLVector tu = camera->u;
    GLVector tv = camera->v;
    glVectorSet(&camera->u, cs*tu.x + sn*tv.x, cs*tu.y + sn*tv.y, cs*tu.z + sn*tv.z);
    glVectorSet(&camera->v, -sn*tu.x + cs*tv.x, -sn*tu.y + cs*tv.y, -sn*tu.z + cs*tv.z);
    glCameraSetCT(camera);
}

void glCameraPitch(GLCamera* camera, float angle)
{

}

void glCameraYaw(GLCamera* camera, float angle)
{
    
}

void glLoadMatrixFromCamera(GLCamera* camera)
{
    glLoadMatrix(camera->CT);
}

void glFrustum(float left, float right, float bottom, float top, float near, float far)
{
    FRUSTUM_M[0] = 2 * near / (right - left);
    FRUSTUM_M[4] = 0;
    FRUSTUM_M[8] = (right + left) / (right - left);
    FRUSTUM_M[12] = 0;
    FRUSTUM_M[1] = 0;
    FRUSTUM_M[5] = 2 * near / (top - bottom);
    FRUSTUM_M[9] = (top + bottom) / (top - bottom);
    FRUSTUM_M[13] = 0;
    FRUSTUM_M[2] = 0;
    FRUSTUM_M[6] = 0;
    FRUSTUM_M[10] = -(far + near) / ((far - near) * 1);
    FRUSTUM_M[14] = -2 * near * far / ((far - near) * 1);
    FRUSTUM_M[3] = 0;
    FRUSTUM_M[7] = 0;
    FRUSTUM_M[11] = -1;
    FRUSTUM_M[15] = 0;
}

void glClear()
{
    memset(IMAGE_DATA, 0xFF000000, sizeof(IMAGE_DATA));
    memset(DEPTHS_DATA, 0xFFFFFFFF, sizeof(DEPTHS_DATA));
}


void newtonC(float points[], float c[], int nd, int count)
{
    float* pmatrix = (float*)malloc(sizeof(float) * (nd - 1) * (1 + count) * count / 2);
    for ( int v = 0; v < (nd - 1); v++ )
    {
        for ( int i = 0; i < count; i++ )
        {
            int n = i - 1;
            int preindex = (count + count + (n - 1) * (-1)) * n / 2;
            n = i;
            int curindex = (count + count + (n - 1) * (-1)) * n / 2;
            for ( int j = i; j < count; j++ )
            {
                int mindex = (curindex + j - i) * (nd - 1) + v;
                NSLog(@"v = %d, mindex = %d", v, mindex);
                if ( i == 0 )
                {
                    pmatrix[mindex] = points[j * nd + v + 1];
                }
                else
                {
                    int index0 = preindex + j - i;
                    int index1 = preindex + j - i + 1;
                    int mindex0 = index0 * (nd - 1) + v;
                    int mindex1 = index1 * (nd - 1) + v;
                    int ptindex0 = (j - i) * nd;
                    int ptindex1 = j * nd;
                    NSLog(@"v = %d, index0 = %d, index1 = %d, pindex0 = %d, pindex1 = %d", v, index0, index1, j-i, j);
                    pmatrix[mindex] = (pmatrix[mindex1] - pmatrix[mindex0]) / (points[ptindex1] - points[ptindex0]);
                }
                NSLog(@"pmatrix[%d] = %f", mindex, pmatrix[mindex]);
                if ( j == i )
                {
                    c[i * (nd - 1) + v] = pmatrix[mindex];
                }
            }
        }
    }
    delete [] pmatrix;
}


void newtonInterpolation(float points[], int nd, int count, float pvs[], int pvcount)
{
    float* c = (float*)malloc(sizeof(float) * count * (nd - 1));
    newtonC(points, c, nd, count);
    for ( int i = 0; i < count; i++ )
    {
        NSLog(@"c%d = %f", i, c[i]);
    }
    float t0 = points[0];
    float t1 = points[(count-1) * nd];
    for ( int v = 0; v < (nd - 1); v++ )
    {
        for ( int i = 0; i < pvcount; i++ )
        {
            int pvindex = i * (nd - 1) + v;
            pvs[pvindex] = 0.0f;
            float t = t0 + (t1 - t0) * i / (pvcount - 1);
            float mt = 1.0f;
            for ( int j = 0; j < count; j++ )
            {
                if ( j > 0 )
                {
                    float tj = points[(j - 1) * nd];
                    mt = mt * (t - tj);
                }
                pvs[pvindex] += c[j * (nd - 1) + v] * mt;
            }
        }
    }
    delete [] c;
}


